#pragma once
/* zarthray - ZarthMath.h 2009
 * The Math used in the Ray Tracing program, many utils live here. */
#include <cstdint>
#include <cstdlib>
#include <cstring>
#include <ostream>
#include <cmath>

using namespace std;

typedef class Vec3f
{
public:
	float	data[4];	//The data stored in the vector, we store 4 floats
						//Because we want to be able to use homogeneous coordinates
public:


	float&	operator[](int position)	//Array operator to make accesing elements
	{									//of the vector easier
		return data[position];
	}

	Vec3f operator *(float f)			//Scalar multiplication
	{
		Vec3f out;
		out[0] = f * data[0];
		out[1] = f * data[1];
		out[2] = f * data[2];
		return out;
	}
	Vec3f operator +(Vec3f rhs)		//addition operator
	{
		Vec3f out;
		out[0] = data[0] + rhs[0];
		out[1] = data[1] + rhs[1];
		out[2] = data[2] + rhs[2];
		return out;
	}

	Vec3f operator -(Vec3f rhs)		//subtraction operator
	{
		return *this + rhs * -1.0f;
	}
	Vec3f(float i = 0.0f, float j = 0.0f, float k = 0.0f)	//default constructor
	{									
		data[0] = i;
		data[1] = j;
		data[2] = k;
		data[3] = 0.0f;
	}

	void normalize()
	{
		float length = sqrtf(data[0]*data[0] + data[1]*data[1] + data[2]*data[2]);
		data[0] = data[0]/length;
		data[1] = data[1]/length;
		data[2] = data[2]/length;
	}

	friend ostream& operator << (ostream& os, const Vec3f v);
	friend float dotProduct(Vec3f a, Vec3f b);
	friend Vec3f crossProduct(Vec3f a, Vec3f b);
	
} Point;


class Vec4f
{
public:
	float data[4];
	float transdata[4];
	Vec4f(float _x, float _y, float _z, float _w)
	{
		data[0] = _x;
		data[1] = _y;
		data[2] = _z;
		data[3] = _w;
	}
	Vec4f()
	{
	}

	operator Vec3f ()
	{
		Vec3f out;
		memcpy(out.data,data,4*sizeof(float));
		out.data[3] = 0;
		return out;
	}
	float& operator[] (unsigned int position)
	{
		return data[position];
	}
	Vec4f(Vec3f v)
	{
		memcpy(data,v.data,4*sizeof(float));
		data[3] = 1.0f;
	}
	Vec4f(float *_data)
	{
		memcpy(data,_data,4*sizeof(float));
	}
	float& operator[](int position)
	{
		return data[position];
	}

};


inline float dotProduct(Vec4f a, Vec4f b)
{
	return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3];
}
inline float dotProduct(Vec3f a, Vec3f b) 
{
	return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}

inline Vec3f crossProduct(Vec3f a, Vec3f b) 
{
	float i = a[1]*b[2] - a[2]*b[1];
	float j = a[2]*b[0] - a[0]*b[2];
	float k = a[0]*b[1] - a[1]*b[0];
	return Vec3f(i,j,k);
}


class Mat4x4f
{
public:
	Vec4f rows[4];
	Vec4f& operator[](int row)
	{
		return rows[row];
	}

	Vec4f col(int _col)
	{
		return Vec4f(rows[0][_col], rows[1][_col], rows[2][_col], rows[3][_col]);
	}

	Mat4x4f& operator *(Mat4x4f& B)
	{
		Mat4x4f out;
		Mat4x4f A = *this;
		for (int i = 0; i < 4; i++)
		{
			for (int j = 0; j < 4; j++)
			{
				out[i][j] = dotProduct(A[i],B.col(j));
			}
		}
	}

	Vec4f& operator *(Vec4f v)
	{
		Vec4f out;
		
		for (int i = 0; i < 4; i++)
		{
			out[i] = dotProduct((*this)[i],v);
		}
	}

};
class uvwBasis
{
	Vec3f u, v, w;
	Mat4x4f T;
	Mat4x4f Tinv;
	Point p;

	uvwBasis(Point _p, Point _at, Vec3f up)
		: p(_p)
	{
		w = _at - p;
		w.normalize();
		u = crossProduct(w,up);
		u.normalize();
		v = crossProduct(u,w);
	}

	//given a vector in xyz, convert to uvw
	Vec3f to_uvw(Vec3f in)
	{
		return T * in;
	}

	//given a vector in uvw, convert to xyz
	Vec3f to_xyz(Vec3f in)
	{
		return Tinv * in;
	}
	uvwBasis()
	{
	}
};


inline ostream& operator << (ostream& os, Vec3f v)	//stream output for vectors
{
	os << "( " << v[0] << ", " << v[1] << ", " << v[2] << ")";
	return os;
}

inline float det(Mat4x4f A)
{
	float a,b,c,d;
	//need the det of the first 3x3 matrix
	/*a = A[0][0]*(A[1][1]*(A[2][2]*A[3][3] - A[2][3]*A[3][2])
				- A[1][2]*(A[2][1]*A[3][3] - A[2][3]*A[3][1])
				+ A[1][3]*(A[2][1]*A[3][2] - A[2][2]*A[3][1]));*/
	float a1 = A[1][1]*(A[2][2]*A[3][3] - A[2][3]*A[3][2]);
	float a2 = A[1][2]*(A[2][1]*A[3][3] - A[2][3]*A[3][1]);
	float a3 = A[1][3]*(A[2][1]*A[3][2] - A[2][2]*A[3][1]);
	a = A[0][0] * (a1 - a2 + a3);
	b = A[0][1]*(A[1][0]*(A[2][2]*A[3][3] - A[2][3]*A[3][2])
				- A[1][2]*(A[2][0]*A[3][3] - A[2][3]*A[3][0])
				+ A[1][3]*(A[2][0]*A[3][2] - A[2][2]*A[3][0]));
	c = A[0][2]*(A[1][0]*(A[2][1]*A[3][3] - A[2][3]*A[3][1])
				- A[1][1]*(A[2][0]*A[3][3] - A[2][3]*A[3][0])
				+ A[1][3]*(A[2][0]*A[3][1] - A[2][1]*A[3][0]));
	d = A[0][3]*(A[1][0]*(A[2][1]*A[3][2] - A[2][2]*A[3][1])
				- A[1][1]*(A[2][0]*A[3][2] - A[2][2]*A[3][0])
				+ A[1][2]*(A[2][0]*A[3][1] - A[2][1]*A[3][0]));

	return a - b + c - d;
}

inline Vec3f operator *(float f, Vec3f v)
{
	return v * f;
}



inline ostream& operator <<(ostream& os, Vec4f v)
{
	os << "[ " << v[0] << ", " << v[1] << ", " << v[2] << ", " << v[3] << "]";
	return os;
}

inline ostream& operator <<(ostream& os, Mat4x4f m)
{
	os << "[" << m[0] << std::endl
		<< m[1] << std::endl
		<< m[2] << std::endl
		<< m[3];
	return os;
}